"""
May 14, 2018

Marine Cadastre Project
www.marinecadastre.gov

NOAA Office for Coastal Management
1234 South Hobson Avenue
Charleston, SC 29405
843-740-1200
ocm.mmc@noaa.gov

Software design and development by:
RPS
55 Village Square Drive
Wakefield, RI 02879
"""

# RK Edits 8/2/2018
# changing line segment algorithm to my modified one
# "C:\Python27\ArcGISx6410.5\python.exe" AIS_TrackBuilder_v3.1_64bit_rk.py "" outworkspace outname

########################################################################################################
#########################################################################################################
#########################################################################################################
## Input Paramaters

#Path to input points (broadcast feature class)
#sInputFile = r"E:\data\temp2.gdb\points_wc2015_01"
#sInputFile = r"E:\data\2014\01_January_2014\Zone8_2014_01.gdb\Broadcast"
sInputFile = r"E:\data\temp2.gdb\test_wc1516_small"

#Path to output workspace (file geodatabase)
sOutWorkspace = r"C:\Users\rklotz\Documents\data\testing\testing1516.gdb"

#name of output trackline feature class
sOutName = "test_small_mp3"

cores = 8 # number of cores to use, 1 sets to serial

### Options for creating tracks ##################################

#ID field name (can keep as is most likely)
sID_Field = "MMSI"

#DateTime field name (can keep as is most likely)
sDT_Field = "BaseDateTime"

#Optional - Attribute fields in the input point data to transfer to the tracklines, the values from the first point for each trackline will be used.
#fields should be separated by a semicolon (;)
#sFieldList = "VesselName;IMO;CallSign;VesselType;Length;Width;Draft;Cargo"
sFieldList = "VesselName;IMO;VesselType;Length;Width;Draft"

# logic for when to break line
##iMaxDist = 20*0.621371       # max distance condition (miles) 
##
### time intervals for other conditions
##iMaxTime1 = 2*60    # connect any within this # of minutes
##iMaxTime2 = 1.9*60  # removing this  # connect if COG close (minutes)
##iMaxTime3 = 24*60   # connect if within predicted cone (minutes)
##
##fSpeedMult = 0.75   # change in speed for prediction of cone
##iCOGrange = 15      # COG range for prediction
##
##iMaxCOGdiff = 15        # difference in COG for close 
##iCOGdiff_predict = 25   # max COG range for matching in predicted cone 
##
##iMaxSOGimplied = 50 # skip any points that implied SOG is very high

##################################################################
#### Unused
#Break Track line method option (comment/uncomment the selected option)
#sLineMethod = "Maximum Time and Distance"
#sLineMethod = "Time Interval"

#Maximum Time in minutes (keep as string in quotes, use 0 if using Time Interval method, default = 30)
#sMaxTimeMinutes = "30"

#Maximum Distance in (statute) miles (keep as string in quotes, use 0 if using Time Interval method, default = 1)
#sMaxDistMiles = "1"

#Time Interval in hours (keep as string in quotes, use 0 if using Maximum Time and Distance method, default = 24)
#sIntervalHours = "0"

#########################################################################################################
#########################################################################################################
#########################################################################################################

import arcpy, os, sys, traceback, datetime, math, re
import multiprocessing
import time
import numpy as np
from functools import partial

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
import warnings


iMaxTimeMinutes = 0
iMaxDistMiles = 0
iIntervalHours = 0
spatialRef_CURSOR = None
lFields = []
#########################################################################################################

class Timer:
    def __init__(self):        
        self.startTime = time.time()
        
    def End(self):
        self.endTime = time.time()
        self.elapTime = self.endTime - self.startTime        
        return self.elapTime

class Progress:
    def __init__(self, total):
        self.iTotal = float(total)
        self.prog = 0
        sys.stdout.write('\r    0%')
        
    def Update(self, current):
        perc = round(current/self.iTotal * 100.0,0)
        if perc > self.prog:
            self.prog = perc
            sys.stdout.write('\r    '+str(perc)+"%")
            sys.stdout.flush()


def FormatCounts(iInCount):
    """Formats various counts to strings with comma separators, for use in dialog messages."""
    sOutCount = re.sub("(\d)(?=(\d{3})+(?!\d))", r"\1,", "%d" % iInCount)
    return sOutCount
 
def CreateDataDict():
    """Process to read through the input point feature class.  Point information is stored in a Python dictionary in memory.
    All point attributes and coordinates (except date/time) are read into Numpy Array.
    The date information is read using the Data Access module search cursor.  
    The two parts are then combined based on the OID into the Python dictionary"""
    
    print "\n  Reading input point data..."
    
    iTotPoints = int(arcpy.GetCount_management(sInputFile).getOutput(0))
    timer = Timer()
    prog = Progress(iTotPoints*3) #once for array, date dict, merge
    tally = 0

    i=0
    fields = ["SHAPE@XY",sID_Field,"OID@"] + ["SOG","COG","Heading"]
    for f in lFields:
        if f.name not in fields:
            fields.append(f.name)
    fields2 = ["OID@",sDT_Field]
        
    if sID_Field == "MMSI":
        where = sID_Field+" > 0"
    else:
        where = None

    #Store all data (except date) into numpy array
    ptData_Array = arcpy.da.FeatureClassToNumPyArray(sInputFile,fields,where,spatialRef_CURSOR,null_value=-9999)
    tally+=iTotPoints
    prog.Update(tally)
     
    #Get date using search cursor, store in temporary dictionary
    dateDict = {}
    dataDicAttrs = {}
    with arcpy.da.SearchCursor(sInputFile,fields2,where,spatialRef_CURSOR,False) as rows:
        for row in rows:      
            sOID = row[0]
            dDateTime = row[1]
            dateDict[sOID] = dDateTime
            
            tally+=1
            prog.Update(tally)
    #to account for points without MMSI
    tally=iTotPoints*2
    prog.Update(tally)
                 
    #combine array and date dictionary into one final dictionary               
    dataDicAIS = {}
    for i in xrange(0,len(ptData_Array[sID_Field])):    
        sID, oid, geo, sog, cog, head = ptData_Array[sID_Field][i], ptData_Array["OID@"][i], ptData_Array["SHAPE@XY"][i] , ptData_Array["SOG"][i], ptData_Array["COG"][i], ptData_Array["Heading"][i]
               
        if dataDicAIS.has_key(sID):
            dataDicAIS[sID][dateDict[oid]] = [geo[0],geo[1], sog, cog, head]
        else:
            dataDicAIS[sID] = {dateDict[oid]:[geo[0],geo[1], sog, cog, head]}
    
        #store the FIRST points attributes for each of the additional fields included  
        #!!! - this data isn't sorted, so the first point may not be the same as the non 64-bit version of this tool.
        if sID not in dataDicAttrs:
            dataDicAttrs[sID] = {}
            for iFieldIndex in range(3,len(fields)):
                val = ptData_Array[fields[iFieldIndex]][i]
                if val == -9999 or val == "-9999":
                    dataDicAttrs[sID][fields[iFieldIndex]] = None
                else: 
                    dataDicAttrs[sID][fields[iFieldIndex]] = val

        tally+=1
        prog.Update(tally)
    
    del sID, oid, geo
    del dateDict, ptData_Array
     
    print "\r    Read "+FormatCounts(iTotPoints)+" points"
    print "    Read Complete:" + str(timer.End())
    del timer, prog
     
    return dataDicAIS, dataDicAttrs, iTotPoints

def SplitDataDictEqually(input_dict, parts, total_points):
    """Splits the point data dictionary into equal parts, for separate processing."""
    return_list = [dict() for idx in xrange(parts)]

    min_pts_per_part = math.ceil(total_points/parts)

    idx = 0
    pts_added = 0
    for k,v in input_dict.iteritems():
        if pts_added < min_pts_per_part:
            return_list[idx][k] = v
            pts_added+=len(v.keys())
        else:
            idx += 1
            pts_added = 0
            
            return_list[idx][k] = v
            pts_added+=len(v.keys())
            
    return return_list

def CreateOuputDataset(parts=0):
    """Creates an empty feature class to store the track polylines.  Adds the required fields, based on the options selected by the user."""
    print "\n  Building output track line feature class..."
     
    tmp = os.path.join('in_memory', 'template')

    has_m = "ENABLED"
    #create a feature class in memory, to use as a template for all feature classes created
    arcpy.CreateFeatureclass_management('in_memory','template',"POLYLINE","",has_m,"DISABLED",sInputFile)

    lOutFieldList = ["SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"]
    arcpy.AddField_management(tmp,sID_Field,"LONG",)
    arcpy.AddField_management(tmp,"TrackStartTime","DATE",)
    arcpy.AddField_management(tmp,"TrackEndTime","DATE",)
     
    #Add other fields
    for f in lFields:
        if f.name not in [sID_Field,sDT_Field]:
            arcpy.AddField_management(tmp,f.name, f.type, f.precision, f.scale, f.length, f.aliasName)
            lOutFieldList.append(f.name)

    #create a temporary FGDB and feature class for each separate process to avoid locking issues
    if parts > 1:
        for i in range(1,parts+1):
            tmp_wkspc_path = os.path.split(sOutWorkspace)[0]
            arcpy.CreateFileGDB_management(tmp_wkspc_path, "temp" + str(i) + ".gdb")
            tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")
            arcpy.CreateFeatureclass_management(tmp_fgdb,sOutName+str(i),"POLYLINE",tmp,has_m,"DISABLED",sInputFile)
    
    #create one feature class since no multiprocessing being used
    else:
        arcpy.CreateFeatureclass_management(sOutWorkspace,sOutName,"POLYLINE",tmp,has_m,"DISABLED",sInputFile)
    
    #delete the temporary template dataset in memory
    arcpy.Delete_management(tmp)
    del tmp
        
    return lOutFieldList

def GetDistance(prevX,prevY,currX,currY):
    """Calculates the distance between two latitude/longitude coordinate pairs, using the Vincenty formula for ellipsoid based distances.
    Returns the distance in miles."""
    #Vincenty Formula
    #Copyright 2002-2012 Chris Veness
    #http://www.movable-type.co.uk/scripts/latlong-vincenty.html
     
    a = 6378137
    b = 6356752.314245
    f = 1/298.257223563

    L = math.radians(currX - prevX)
    U1 = math.atan((1-f)*math.tan(math.radians(prevY)))
    U2 = math.atan((1-f)*math.tan(math.radians(currY)))
    sinU1 = math.sin(U1)
    sinU2 = math.sin(U2)
    cosU1 = math.cos(U1)
    cosU2 = math.cos(U2)

    lam = L
    lamP = 9999999999
    iter_count = 0
    while abs(lam-lamP) > 1e-12 and iter_count <= 100:
        sinLam = math.sin(lam)
        cosLam = math.cos(lam)
        sinSigma = math.sqrt((cosU2*sinLam)*(cosU2*sinLam) + (cosU1*sinU2-sinU1*cosU2*cosLam)*(cosU1*sinU2-sinU1*cosU2*cosLam))

        if sinSigma == 0:
            return 0

        cosSigma = sinU1*sinU2+cosU1*cosU2*cosLam
        sigma = math.atan2(sinSigma,cosSigma)
          
        sinAlpha = cosU1*cosU2*sinLam/sinSigma
        cosSqAlpha = 1 - sinAlpha*sinAlpha
        if cosSqAlpha == 0: #catches zero division error if points on equator
            cos2SigmaM = 0
        else:
            cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha

        if cos2SigmaM == None:
            cos2SigmaM = 0

        C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        lamP = lam
        lam = L+(1-C)*f*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))

        iter_count+=1

    uSq = cosSqAlpha*(a*a - b*b)/(b*b)
    A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B = uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)))
    deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
    s = b*A*(sigma-deltaSigma)

    #convert s in meters to miles
    s_miles = s*0.0006213712
    return s_miles

def GetTimeDifference(prevDT,currDT):
    """Calculates the difference between to date and time variables. The difference is returned in minutes."""
    timeDelta = currDT - prevDT
    #totMinutes = (timeDelta.days*1440) + (timeDelta.seconds/60)
    totMinutes = timeDelta.total_seconds() / 60 # a better option

    return totMinutes

def MultiProcess_Main(dataDicAIS, dataDicAttrs, lOutFieldList, iTotPoints,cores):
    """Main process to split the data and run seperate processes to build the track lines for each available CPU."""
    print "\n  Generating track lines (Multiprocess)..."
    
    timer = Timer()        
    #parts = multiprocessing.cpu_count() - 1
    parts = cores   
    prog = Progress(1 + parts + 2) #split, each sub-process, merge, delete temp
    tally = 0
     
    tmp_wkspc_path = os.path.split(sOutWorkspace)[0]
    
    #split the point data into separate dictionaries for each process.
    split_dicts = SplitDataDictEqually(dataDicAIS, parts, iTotPoints)
    tally+=1
    prog.Update(tally)
    
    #start each separate process asynchronously, and wait for all to finish before moving on.
    pool = multiprocessing.Pool(parts)
    jobs = []
    for i in range(1,parts+1):
        tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")
        sTempFC = tmp_fgdb+"\\"+sOutName+str(i)
        splitDict = split_dicts[i-1]
        
        tally+=1
        p_callback = partial(MultiProcess_Status,p=prog,t=tally)
     
        #if sLineMethod == "Maximum Time and Distance":
        jobs.append(pool.apply_async(AddLines_Segmented, (splitDict,dataDicAttrs,lOutFieldList,sTempFC,spatialRef_CURSOR), callback=p_callback))
        #jobs.append(pool.apply_async(AddLines_Segmented, (splitDict,dataDicAttrs,lOutFieldList,sTempFC,spatialRef_CURSOR) ))
        #elif sLineMethod == "Time Interval":
            #jobs.append(pool.apply_async(AddLines_Interval, (splitDict,dataDicAttrs,lOutFieldList,sTempFC,iIntervalHours,spatialRef_CURSOR), callback=p_callback))

    pool.close()
    pool.join()

    del split_dicts

    #Merge the temporary output feature classes in separate FGDBs into one feature class in the selected output geodatabase
    to_merge = []
    for i in range(1,parts+1):
        tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")
        to_merge.append(tmp_fgdb+"\\"+sOutName+str(i))

    arcpy.Merge_management(to_merge, sOutWorkspace+"\\"+sOutName)
    tally+=1
    prog.Update(tally)

    #Delete the temporary geodatabases
    for i in range(1,parts+1):
        tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")
        arcpy.Delete_management(tmp_fgdb)
    tally+=1
    prog.Update(tally)
        
    iTotTracks = int(arcpy.GetCount_management(sOutWorkspace+"\\"+sOutName).getOutput(0))
    print "\r    Unique "+sID_Field+"'s= "+FormatCounts(len(dataDicAIS))
    print "    Track lines created = "+FormatCounts(iTotTracks)
    print "    Track lines created in:" + str(timer.End())
    
    del timer
    
def MultiProcess_Status(x, p, t):
    """Simple function to get the status of each separate track line building process."""
    p.Update(t)    

def ProjectPoint(pnt,c,s,t) :
    # calculates prediction of next point based on cog (c) and speed(s) and time gap (t)
    # c is angle from north, s is knots; t is in hours
    # returns values in decimal degrees

    # convert c (cog) to angle from x-axis
    c90 = 90 - c  ;
    if c90<0 :
        c90 = 360 + c90  

    a = c90*math.pi / 180 ; # convert to radians

    #print a
    #print t*s*math.cos(a)
    #print t*s*math.sin(a)

    #pnt_prj  = pnt
    #print pnt_prj

    pnt_prj = [0]*2
    # dividing by 60 to convert from nm to degrees
    pnt_prj[0] = pnt[0] + ( ( t*s*math.cos(a) ) / 60 ) 
    pnt_prj[1] = pnt[1] + ( ( t*s*math.sin(a) ) / 60 ) 

    #print pnt_prj

    return pnt_prj

def degDiff(np_arr) :
    # calculating difference between two angles
    # catching and ignoring warnings
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        normDeg = np.mod( ( np_arr[1] - np_arr[0]) , 360 )
        deg_diff = min(360-normDeg, normDeg)

    return deg_diff
 
def AddLines_Segmented(dataDicAIS,dataDicAttrs,lOutFieldList,sOutputFC,spatialRef_CURSOR,progress=False):
    """Creates track polylines and saves them in the ouput feature class.  Track lines are created using a maximum distance and
    maximum time threshold between points."""

    # logic for when to break line
    iMaxDist = 20*0.621371       # max distance condition (miles) 

    # time intervals for other conditions
    iMaxTime1 = 2*60    # connect any within this # of minutes
    iMaxTime2 = 1.9*60  # removing this  # connect if COG close (minutes)
    iMaxTime3 = 24*60   # connect if within predicted cone (minutes)

    fSpeedMult = 0.75   # change in speed for prediction of cone
    iCOGrange = 15      # COG range for prediction

    iMaxCOGdiff = 15        # difference in COG for close 
    iCOGdiff_predict = 25   # max COG range for matching in predicted cone 

    iMaxSOGimplied = 50 # skip any points that implied SOG is very high

    
    prog = None
    tally = 0
    if progress:
        prog = Progress(len(dataDicAIS))

    #cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))
    cur = arcpy.da.InsertCursor(sOutputFC, lOutFieldList)

    pnt1 = arcpy.Point()
    pnt2 = arcpy.Point()
    lineArray = arcpy.Array()

    for key in sorted(dataDicAIS.iterkeys()):
        tally+=1
        iSegCount = 0

        lAddAttrs = []
        for f in lOutFieldList:
            if f not in ["SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"]:
                if dataDicAttrs[key].has_key(f):
                    lAddAttrs.append(dataDicAttrs[key][f])
                else:
                    lAddAttrs.append(None)

        #list of date time dictionary keys
        dtTimes = sorted(dataDicAIS[key].keys())

        #start and end date variables: start date set with the start of any new line.  end date constantly updated with the second point
        #so that whenever the line ends the last datetime will already be set.
        dtSegStart = None
        dtSegEnd = None

        ### RK CHANGES 
        # get list of SOG and COG
        SOGall = [ dataDicAIS[key][k][2] for k in dtTimes ]
        COGall = [ dataDicAIS[key][k][3] for k in dtTimes ]

        # checking if SOG and COG are misreported (multiplied by 10)
        # SOG too high if any value greater than XX
        # COG too high if any value greater than 360 
        SOGhigh = any(i > 1020 for i in SOGall)
        COGhigh = any(i > 360 for i in COGall)

        # if either one is high then divide both COG and SOG by 10
        valMult = 1 
        if SOGhigh or COGhigh :
            valMult = 1/10

        #skip MMSI if there is only one point provided            
        if len(dtTimes) > 1:
            #for i in range(0,len(dtTimes)-1):
            i1 = 0# index of first point
            for i2 in range(1,len(dtTimes)):
                                
                # getting SOG and COG as numpy arrays
                SOGvec = np.array([dataDicAIS[key][dtTimes[i1]][2]*valMult , dataDicAIS[key][dtTimes[i2]][2]*valMult],dtype=np.dtype('f8'))
                COGvec = np.array([dataDicAIS[key][dtTimes[i1]][3]*valMult , dataDicAIS[key][dtTimes[i2]][3]*valMult],dtype=np.dtype('f8'))
                HEADvec = np.array([ dataDicAIS[key][dtTimes[i1]][4] , dataDicAIS[key][dtTimes[i2]][4] ],dtype=np.dtype('f8'))

                # replacing with NaN if bad
                SOGvec[SOGvec == 102] = np.nan
                COGvec[COGvec == 360] = np.nan
                HEADvec[COGvec == 511] = np.nan
                
                pnt1.X = dataDicAIS[key][dtTimes[i1]][0]
                pnt1.Y = dataDicAIS[key][dtTimes[i1]][1]
                pnt1.M = time.mktime( dtTimes[i1].timetuple() )
                
                pnt2.X = dataDicAIS[key][dtTimes[i2]][0]
                pnt2.Y = dataDicAIS[key][dtTimes[i2]][1]
                pnt2.M = time.mktime( dtTimes[i2].timetuple() )
                
                distance = GetDistance(pnt1.X,pnt1.Y,pnt2.X,pnt2.Y)
                timeDiff = GetTimeDifference(dtTimes[i1],dtTimes[i2])

                SOG_implied = ( distance / (timeDiff/60) ) * 0.868976  # implied speed over ground (mph to knots)
   
                # getting average SOG and implied COG
                # ignoring warning if both nans
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore", category=RuntimeWarning)
                    SOG_reported_avg = np.nanmean(SOGvec)   # avg reported SOG

                angle_implied = math.atan2( pnt2.Y - pnt1.Y , pnt2.X - pnt1.X )  * 180 / math.pi      
                COG_implied = angle_implied + 90   ;
                if COG_implied>360 :
                    COG_implied = COG_implied - 360   

                # absolute value of COG difference and Heading difference
                cogDiff = degDiff(COGvec)
                headDiff = degDiff(HEADvec)

                # break line unless one of the following is true
                continueFlag = False

                if distance<=iMaxDist :
                    continueFlag = True
                
                elif timeDiff <= iMaxTime1 :
                    continueFlag = True

                elif timeDiff <= iMaxTime2 :
                    # if either cogDiff or headDiff is within bound
                    if cogDiff<=iMaxCOGdiff or headDiff<=iMaxCOGdiff :
                        continueFlag = True
                        
                elif timeDiff <= iMaxTime3 :
                    # calculate projected points (if a valid COG or heading)

                    # get course or heading
                    c = COGvec[0]
                    if np.isnan(c) :
                        c = HEADvec[0]

                    # only check if valid c
                    if not np.isnan(c) :
                        pnt_start = [ pnt1.X , pnt1.Y ]    
                        s = SOG_implied
                        t = timeDiff/60
                        
                        c_minus = c - iCOGrange
                        if c_minus<0 :
                            c_minus = 360 + c_minus

                        c_plus = c + iCOGrange
                        if c_plus>360:
                            c_plus = c_plus - 360

                        # Making cone
                        # speed*(1+y)
                        pnt_pred_1 = ProjectPoint(pnt_start,c_plus,s*(1+fSpeedMult),t) # cog + x 
                        pnt_pred_2 = ProjectPoint(pnt_start,c,s*(1+fSpeedMult),t) 
                        pnt_pred_3 = ProjectPoint(pnt_start,c_minus,s*(1+fSpeedMult),t) # cog - x

                        point2 = Point(pnt2.X,pnt2.Y)
                        poly_range = Polygon( [ pnt_start, pnt_pred_1, pnt_pred_2, pnt_pred_3 ] )

                        # continue if within range of predicted point
                        if poly_range.contains(point2) :
                            if cogDiff<=iCOGdiff_predict or headDiff<=iCOGdiff_predict : 
                                continueFlag = True
                
                # skipping major anomolies (very fast implied SOG)
                skipFlag = False
                if SOG_implied > iMaxSOGimplied : 
                    skipFlag=True
                    
                # update line   
                if continueFlag :
                    if not skipFlag :
                        if iSegCount == 0:
                            dtSegStart = dtTimes[i1]
                            lineArray.add(pnt1)
                            lineArray.add(pnt2)

                        else:
                            lineArray.add(pnt2)

                        dtSegEnd = dtTimes[i2]
                        iSegCount+=1

                #Break the line when continue flag is false
                else:                   
                    if lineArray.count > 1:
                        polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
                        cur.insertRow((polyline,key,dtSegStart,dtSegEnd) + tuple(lAddAttrs))
                        del polyline
           
                    lineArray.removeAll()
                    dtSegStart = None
                    dtSegEnd = None

                    ##reset seg count,since line ended early due to thresholds
                    iSegCount = 0

                # iterate i1 to i2 (unless point is skipped)
                if not skipFlag :
                    i1=i2 

            #end line since end of this MMSI set
            if lineArray.count > 1:
                polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
                cur.insertRow((polyline,key,dtSegStart,dtSegEnd) + tuple(lAddAttrs))
                del polyline

            lineArray.removeAll()
            dtSegStart = None
            dtSegEnd = None
            ########
        
        if progress:
            prog.Update(tally)

    del prog
    del cur

def AddLines_Interval(dataDicAIS,dataDicAttrs,lOutFieldList,sOutputFC,iIntervalHours,spatialRef_CURSOR,progress=False):
    """Creates track polylines and saves them in the ouput feature class.  Track lines are created based on the specified time interval."""
    
    prog = None
    tally = 0
    if progress:
        prog = Progress(len(dataDicAIS))
        
    #cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))
    cur = arcpy.da.InsertCursor(sOutputFC, lOutFieldList)

    pnt1 = arcpy.Point()
    pnt2 = arcpy.Point()
    lineArray = arcpy.Array()

    #create time interval from input parameter hours
    tdInterval = datetime.timedelta(hours=iIntervalHours)

    for key in sorted(dataDicAIS.iterkeys()):
        tally+=1
        iSegCount = 0
        
        lAddAttrs = []
        for f in lOutFieldList:
            if f not in ["SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"]:
                if dataDicAttrs[key].has_key(f):
                    lAddAttrs.append(dataDicAttrs[key][f])
                else:
                    lAddAttrs.append(None)

        #list of date time dictionary keys
        dtTimes = sorted(dataDicAIS[key].keys())

        #start and end date variables: start date set with the start of any new line.  end date constatnly updated with the second point
        #so that whenever the line ends the last datetime will already be set.
        dtSegStart = None
        dtSegEnd = None

        #initialize time interval using the first point as the start, and the end based on the selected interval
        #updated as the interval threshold is reached.
        dtIntervalBegin = dtTimes[0]
        dtIntervalEnd = dtIntervalBegin + tdInterval

        #skip MMSI if there is only one point provided
        if len(dtTimes) > 1:
            for i in range(0,len(dtTimes)-1):                    
                pnt1.X = dataDicAIS[key][dtTimes[i]][0]
                pnt1.Y = dataDicAIS[key][dtTimes[i]][1]

                pnt2.X = dataDicAIS[key][dtTimes[i+1]][0]
                pnt2.Y = dataDicAIS[key][dtTimes[i+1]][1]

                if dtTimes[i] >= dtIntervalBegin and dtTimes[i+1] <= dtIntervalEnd:
                    if iSegCount == 0:
                        dtSegStart = dtTimes[i]
                        lineArray.add(pnt1)
                        lineArray.add(pnt2)
                    else:
                        lineArray.add(pnt2)

                    dtSegEnd = dtTimes[i+1]
                    iSegCount+=1

                #Break the line as the next point is outside the defined interval
                else:
                    if lineArray.count > 1:                              
                        polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)
                        cur.insertRow((polyline,key,dtSegStart,dtSegEnd) + tuple(lAddAttrs))
                        del polyline
           
                    lineArray.removeAll()
                    dtSegStart = None
                    dtSegEnd = None

                    #update for next time interval
                    dtIntervalBegin = dtIntervalEnd
                    dtIntervalEnd = dtIntervalBegin + tdInterval

                    ##reset seg count,since line ended early due to thresholds
                    iSegCount = 0

            #end line since end of this MMSI set
            if lineArray.count > 1:
                polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)
                cur.insertRow((polyline,key,dtSegStart,dtSegEnd) + tuple(lAddAttrs))
                del polyline
                    
            lineArray.removeAll()               
            dtSegStart = None
            dtSegEnd = None
            dtIntervalBegin = None
            dtIntervalEnd = None

        if progress:
            prog.Update(tally)

    del prog
    del cur

def CreateTracks():
    """Main function that calls all of the other functions.  Determines which functions are run, based on the user selected parameters."""
    #Create Dictionary of broadcast data Key = ID, Value = Dictionary of {DateTime:[x,y,voyageID]}
    dataDicAIS, dataDicAttrs, iTotPoints = CreateDataDict()

    sOutputFC = sOutWorkspace+"\\"+sOutName
     
    #If there are more than 2 CPUs then use multiprocessing approach, otherwise only use 1
    #  leaving one CPU to run operating system.
    #if multiprocessing.cpu_count() <= 2:
    if cores<2 :
        #Create output feature class          
        lOutFieldList = CreateOuputDataset(0)
            
        #build track lines
        #print "\n  Generating track lines..."
        #if sLineMethod == "Maximum Time and Distance":     
        AddLines_Segmented(dataDicAIS,dataDicAttrs,lOutFieldList,sOutputFC,spatialRef_CURSOR,True)        
        #elif sLineMethod == "Time Interval":
        #    AddLines_Interval(dataDicAIS,dataDicAttrs,lOutFieldList,sOutputFC,iIntervalHours,spatialRef_CURSOR,True)            
            
        iTotTracks = int(arcpy.GetCount_management(sOutputFC).getOutput(0))
        print "\r    Unique "+sID_Field+"'s= "+FormatCounts(len(dataDicAIS))
        print "    Track lines created = "+FormatCounts(iTotTracks)
     
    else:
        #Create multiple temproary output feature classes
        #lOutFieldList = CreateOuputDataset(multiprocessing.cpu_count() - 1)
        lOutFieldList = CreateOuputDataset(cores)
          
        #Use multiprocessing method to create track lines 
        MultiProcess_Main(dataDicAIS,dataDicAttrs,lOutFieldList,iTotPoints,cores)
         
    print "\n"

def CheckParameters():
    """This function performs the checks required to make sure all the input parameters were provided and accurate.
    These are the same checks performed within the Toolbox script tool validation code."""
    global spatialRef_CURSOR
    global iMaxTimeMinutes, iMaxDistMiles, iIntervalHours
      
    print "\n  Checking Input Parameters..."

    
    bPass = True
    Messages = []
    bRealAISdb = False
    desc = None
    
    #input points checks
    if arcpy.Exists(sInputFile) == False:
        bPass = False
        Messages.append("Input Points do not exist")
    else:
        desc = arcpy.Describe(sInputFile)
        if desc.datasetType <> "FeatureClass":
            bPass = False
            Messages.append("Input is not a feature class")
        else:
            if desc.shapeType <> "Point":
                bPass = False
                Messages.append("Input is not a point feature class")
            
            else:
                prjFile = os.path.join(arcpy.GetInstallInfo()["InstallDir"],"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj")
                spatialRef_WGS84 = arcpy.SpatialReference(prjFile)
                spatialRef_input = arcpy.Describe(sInputFile).spatialReference

                if spatialRef_WGS84.name <> spatialRef_input.name:
                    spatialRef_CURSOR = spatialRef_WGS84
                else:
                    spatialRef_CURSOR = None
        
                fields = desc.fields
                #ID Field checks
                bFieldExists = False
                for field in fields:
                    if field.name == sID_Field:
                        bFieldExists = True
                        if field.type not in ["SmallInteger","Integer"]:
                            bPass = False
                            Messages.append("Input ID field is not an integer field")
                        break
                if bFieldExists == False:
                    bPass = False
                    Messages.append("Input ID field does not exist in input feature class")                        
            
                #DateTime Field checks
                bFieldExists = False
                for field in fields:
                    if field.name == sDT_Field:
                        bFieldExists = True
                        if field.type <> "Date":
                            bPass = False
                            Messages.append("Input Date/Time field is not an date field")
                        break
                if bFieldExists == False:
                    bPass = False
                    Messages.append("Input Date/Time field does not exist in input feature class")                
                    
                #Check list of fields (type, exists, other?)
                if sFieldList <> "":
                    try:
                        lFieldNames = sFieldList.split(";")
                        for sFName in lFieldNames:
                            bFieldExists = False
                            for field in fields:
                                if field.name == sFName:
                                    bFieldExists = True
                                    if field.type not in ['SmallInteger', 'Integer', 'Single', 'Double', 'String']:
                                        bPass = False
                                        Messages.append("Input field ("+sFName+") is not a valid field type ('SmallInteger', 'Integer', 'Single', 'Double', 'String')")
                                    else:
                                        lFields.append(field)
                                    break
                            if bFieldExists == False:
                                bPass = False
                                Messages.append("Input field ("+sFName+") does not exist in input feature class")
                            
                    except:
                        bPass = False
                        Messages.append("Input list of fields to transfer to the tracklines is invalid")                    

    #Output Workspace checks
    if arcpy.Exists(sOutWorkspace) == False:
        bPass = False
        Messages.append("Output workspace does not exist")
    else:
        desc = arcpy.Describe(sOutWorkspace)
        if desc.workspaceType <> "LocalDatabase":
            bPass = False
            Messages.append("Output workspace is not a geodatabase")

    #Check output name
    if sOutName <> str(arcpy.ValidateTableName(sOutName,sOutWorkspace)):
        bPass = False
        Messages.append("Output trackline feature class name is invalid")
    
    #Check Break trackline option
    #if sLineMethod not in ["Maximum Time and Distance","Time Interval"]:
    #    bPass = False
    #    Messages.append("Break trackline method parameter is not a valid option")
        
##    #Check Maximum Time in minutes
##    try:
##        if sLineMethod == "Maximum Time and Distance":            
##            test = int(sMaxTimeMinutes)
##            if int(sMaxTimeMinutes) <= 0:
##                bPass = False
##                Messages.append("Maximum Time must be greater than zero when using the Maximum Time and Distance option")
##            else:
##                iMaxTimeMinutes = int(sMaxTimeMinutes)
##    except:
##        bPass = False
##        Messages.append("Maximum Time is not numeric")
##                
##    #Check Maximum Distance
##    try:
##        if sLineMethod == "Maximum Time and Distance":
##            test = int(sMaxDistMiles)
##            if int(sMaxDistMiles) <= 0:
##                bPass = False
##                Messages.append("Maximum Distance must be greater than zero when using the Maximum Time and Distance option")
##            else:
##                iMaxDistMiles = int(sMaxDistMiles)
##    except:
##        bPass = False
##        Messages.append("Maximum Distance is not numeric")    
##
##    #Check Time Interval
##    try:
##        if sLineMethod == "Time Interval":            
##            test = int(sIntervalHours)
##            if int(sIntervalHours) <= 0:
##                bPass = False
##                Messages.append("Time Interval must be greater than zero when using the Time Interval option")
##            else:
##                iIntervalHours = int(sIntervalHours)
##    except:
##        bPass = False
##        Messages.append("Time Interval is not numeric")   
        
    return bPass, Messages

def HelpText(Messages = []):
    """Returns the command line syntax and examples of usage, when there was a mistake."""
    if len(Messages) > 0:
        print "\n\n"
        for message in Messages:
            print "  ERROR: " + message        
    
    if len(sys.argv) <> 1:
        print "\n\n"    
        print "Trackbuilder Command Line Syntax:"
        print """<path to x64 python.exe> <path to Trackbuilder .py script> <Path to input points> <Path to output workspace> <name of output trackline feature class>"""
    
        print ""
        print "Examples:"
        print """"C:\Python27\ArcGISx6410.6\python.exe" "C:\AIS\AIS_TrackBuilder_v3.1_64bit.py" "C:\AIS\Zone19_2011_07.gdb\Broadcast" "C:\AIS\Zone19_2011_07.gdb" "TrackLine" """
        
    
##########################################
##########################################
if __name__ == '__main__':
    
    timer = Timer()

    bRun = True
    Errors = []
   
    dicInstallInfo = arcpy.GetInstallInfo()
    fArcVersion = float(dicInstallInfo["Version"].split(".")[0] + "." + dicInstallInfo["Version"].split(".")[1])
    if fArcVersion > 10.1:            
        if sys.exec_prefix.find("x64") > 0:
            if len(sys.argv) == 1:
                bRun, Errors = CheckParameters()    
            else:
                #only proceed if all parameters are present
                if len(sys.argv) <> 4 :
                    bRun = False
                else:
                    sInputFile = sys.argv[1]
                    #sID_Field = sys.argv[2]
                    #sDT_Field = sys.argv[3]
                    sOutWorkspace = sys.argv[2]
                    sOutName = sys.argv[3]
                    #sLineMethod = sys.argv[6]
                    #sMaxTimeMinutes = sys.argv[7]
                    #sMaxDistMiles = sys.argv[8]
                    #sIntervalHours = sys.argv[9]
                    
                    #if len(sys.argv) == 11:
                    #    sFieldList = sys.argv[10] 
                    #else:
                    #    sFieldList = ""

                    print sInputFile, sOutWorkspace, sOutName
                                
                    bRun, Errors = CheckParameters()
        else:
            bRun = False
            msg = "   ERROR: The x64 version of the TrackBuilder can only be used with a 64bit version of Python"
            print msg
    else:
        bRun = False
        msg = "   ERROR: TrackBuilder can only be used with ArcGIS Desktop 10.1 or later"
        print msg
        

    
    #Run the process if all parameters were validated
    if bRun:
        try:
            arcpy.env.overwriteOutput = True

            CreateTracks()
            
            print "Track Lines complete!"
            print "  Process complete in " + str(datetime.timedelta(seconds=timer.End()))
            del timer
             
        except arcpy.ExecuteError:
            for msg in range(0, arcpy.GetMessageCount()):
                if arcpy.GetSeverity(msg) == 2:
                    print arcpy.GetMessage(msg)
         
        except:
            tb = sys.exc_info()[2]
            for i in range(0,len(traceback.format_tb(tb))):
                tbinfo = traceback.format_tb(tb)[i]
                msg = "   ERROR: \n   " + tbinfo + "   Error Info:\n   " + str(sys.exc_info()[1])

            print msg

        finally :
            print "Done"
    
    #show syntax and examples if parameters were missing or invalid.
    else:
        HelpText(Errors)
